home *** CD-ROM | disk | FTP | other *** search
/ IRIX 6.2 Development Libraries / SGI IRIX 6.2 Development Libraries.iso / dist / complib.idb / usr / share / catman / p_man / cat3 / complib / clatrs.z / clatrs
Text File  |  1996-03-14  |  8KB  |  265 lines
  1.  
  2.  
  3.  
  4. CCCCLLLLAAAATTTTRRRRSSSS((((3333FFFF))))                                                          CCCCLLLLAAAATTTTRRRRSSSS((((3333FFFF))))
  5.  
  6.  
  7.  
  8. NNNNAAAAMMMMEEEE
  9.      CLATRS - solve one of the triangular systems   A * x = s*b, A**T * x =
  10.      s*b, or A**H * x = s*b,
  11.  
  12. SSSSYYYYNNNNOOOOPPPPSSSSIIIISSSS
  13.      SUBROUTINE CLATRS( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE, CNORM,
  14.                         INFO )
  15.  
  16.          CHARACTER      DIAG, NORMIN, TRANS, UPLO
  17.  
  18.          INTEGER        INFO, LDA, N
  19.  
  20.          REAL           SCALE
  21.  
  22.          REAL           CNORM( * )
  23.  
  24.          COMPLEX        A( LDA, * ), X( * )
  25.  
  26. PPPPUUUURRRRPPPPOOOOSSSSEEEE
  27.      CLATRS solves one of the triangular systems
  28.  
  29.      with scaling to prevent overflow.  Here A is an upper or lower triangular
  30.      matrix, A**T denotes the transpose of A, A**H denotes the conjugate
  31.      transpose of A, x and b are n-element vectors, and s is a scaling factor,
  32.      usually less than or equal to 1, chosen so that the components of x will
  33.      be less than the overflow threshold.  If the unscaled problem will not
  34.      cause overflow, the Level 2 BLAS routine CTRSV is called. If the matrix A
  35.      is singular (A(j,j) = 0 for some j), then s is set to 0 and a non-trivial
  36.      solution to A*x = 0 is returned.
  37.  
  38.  
  39. AAAARRRRGGGGUUUUMMMMEEEENNNNTTTTSSSS
  40.      UPLO    (input) CHARACTER*1
  41.              Specifies whether the matrix A is upper or lower triangular.  =
  42.              'U':  Upper triangular
  43.              = 'L':  Lower triangular
  44.  
  45.      TRANS   (input) CHARACTER*1
  46.              Specifies the operation applied to A.  = 'N':  Solve A * x = s*b
  47.              (No transpose)
  48.              = 'T':  Solve A**T * x = s*b  (Transpose)
  49.              = 'C':  Solve A**H * x = s*b  (Conjugate transpose)
  50.  
  51.      DIAG    (input) CHARACTER*1
  52.              Specifies whether or not the matrix A is unit triangular.  = 'N':
  53.              Non-unit triangular
  54.              = 'U':  Unit triangular
  55.  
  56.      NORMIN  (input) CHARACTER*1
  57.              Specifies whether CNORM has been set or not.  = 'Y':  CNORM
  58.              contains the column norms on entry
  59.              = 'N':  CNORM is not set on entry.  On exit, the norms will be
  60.  
  61.  
  62.  
  63.                                                                         PPPPaaaaggggeeee 1111
  64.  
  65.  
  66.  
  67.  
  68.  
  69.  
  70. CCCCLLLLAAAATTTTRRRRSSSS((((3333FFFF))))                                                          CCCCLLLLAAAATTTTRRRRSSSS((((3333FFFF))))
  71.  
  72.  
  73.  
  74.              computed and stored in CNORM.
  75.  
  76.      N       (input) INTEGER
  77.              The order of the matrix A.  N >= 0.
  78.  
  79.      A       (input) COMPLEX array, dimension (LDA,N)
  80.              The triangular matrix A.  If UPLO = 'U', the leading n by n upper
  81.              triangular part of the array A contains the upper triangular
  82.              matrix, and the strictly lower triangular part of A is not
  83.              referenced.  If UPLO = 'L', the leading n by n lower triangular
  84.              part of the array A contains the lower triangular matrix, and the
  85.              strictly upper triangular part of A is not referenced.  If DIAG =
  86.              'U', the diagonal elements of A are also not referenced and are
  87.              assumed to be 1.
  88.  
  89.      LDA     (input) INTEGER
  90.              The leading dimension of the array A.  LDA >= max (1,N).
  91.  
  92.      X       (input/output) COMPLEX array, dimension (N)
  93.              On entry, the right hand side b of the triangular system.  On
  94.              exit, X is overwritten by the solution vector x.
  95.  
  96.      SCALE   (output) REAL
  97.              The scaling factor s for the triangular system A * x = s*b,  A**T
  98.              * x = s*b,  or  A**H * x = s*b.  If SCALE = 0, the matrix A is
  99.              singular or badly scaled, and the vector x is an exact or
  100.              approximate solution to A*x = 0.
  101.  
  102.      CNORM   (input or output) REAL array, dimension (N)
  103.  
  104.              If NORMIN = 'Y', CNORM is an input argument and CNORM(j) contains
  105.              the norm of the off-diagonal part of the j-th column of A.  If
  106.              TRANS = 'N', CNORM(j) must be greater than or equal to the
  107.              infinity-norm, and if TRANS = 'T' or 'C', CNORM(j) must be
  108.              greater than or equal to the 1-norm.
  109.  
  110.              If NORMIN = 'N', CNORM is an output argument and CNORM(j) returns
  111.              the 1-norm of the offdiagonal part of the j-th column of A.
  112.  
  113.      INFO    (output) INTEGER
  114.              = 0:  successful exit
  115.              < 0:  if INFO = -k, the k-th argument had an illegal value
  116.  
  117. FFFFUUUURRRRTTTTHHHHEEEERRRR DDDDEEEETTTTAAAAIIIILLLLSSSS
  118.      A rough bound on x is computed; if that is less than overflow, CTRSV is
  119.      called, otherwise, specific code is used which checks for possible
  120.      overflow or divide-by-zero at every operation.
  121.  
  122.      A columnwise scheme is used for solving A*x = b.  The basic algorithm if
  123.      A is lower triangular is
  124.  
  125.           x[1:n] := b[1:n]
  126.  
  127.  
  128.  
  129.                                                                         PPPPaaaaggggeeee 2222
  130.  
  131.  
  132.  
  133.  
  134.  
  135.  
  136. CCCCLLLLAAAATTTTRRRRSSSS((((3333FFFF))))                                                          CCCCLLLLAAAATTTTRRRRSSSS((((3333FFFF))))
  137.  
  138.  
  139.  
  140.           for j = 1, ..., n
  141.                x(j) := x(j) / A(j,j)
  142.                x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j]
  143.           end
  144.  
  145.      Define bounds on the components of x after j iterations of the loop:
  146.         M(j) = bound on x[1:j]
  147.         G(j) = bound on x[j+1:n]
  148.      Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}.
  149.  
  150.      Then for iteration j+1 we have
  151.         M(j+1) <= G(j) / | A(j+1,j+1) |
  152.         G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] |
  153.                <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | )
  154.  
  155.      where CNORM(j+1) is greater than or equal to the infinity-norm of column
  156.      j+1 of A, not counting the diagonal.  Hence
  157.  
  158.         G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | )
  159.                      1<=i<=j
  160.      and
  161.  
  162.         |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| )
  163.                                       1<=i< j
  164.  
  165.      Since |x(j)| <= M(j), we use the Level 2 BLAS routine CTRSV if the
  166.      reciprocal of the largest M(j), j=1,..,n, is larger than
  167.      max(underflow, 1/overflow).
  168.  
  169.      The bound on x(j) is also used to determine when a step in the columnwise
  170.      method can be performed without fear of overflow.  If the computed bound
  171.      is greater than a large constant, x is scaled to prevent overflow, but if
  172.      the bound overflows, x is set to 0, x(j) to 1, and scale to 0, and a
  173.      non-trivial solution to A*x = 0 is found.
  174.  
  175.      Similarly, a row-wise scheme is used to solve A**T *x = b  or A**H *x =
  176.      b.  The basic algorithm for A upper triangular is
  177.  
  178.           for j = 1, ..., n
  179.                x(j) := ( b(j) - A[1:j-1,j]' * x[1:j-1] ) / A(j,j)
  180.           end
  181.  
  182.      We simultaneously compute two bounds
  183.           G(j) = bound on ( b(i) - A[1:i-1,i]' * x[1:i-1] ), 1<=i<=j
  184.           M(j) = bound on x(i), 1<=i<=j
  185.  
  186.      The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we add
  187.      the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1.  Then the
  188.      bound on x(j) is
  189.  
  190.           M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) |
  191.  
  192.  
  193.  
  194.  
  195.                                                                         PPPPaaaaggggeeee 3333
  196.  
  197.  
  198.  
  199.  
  200.  
  201.  
  202. CCCCLLLLAAAATTTTRRRRSSSS((((3333FFFF))))                                                          CCCCLLLLAAAATTTTRRRRSSSS((((3333FFFF))))
  203.  
  204.  
  205.  
  206.                <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| )
  207.                          1<=i<=j
  208.  
  209.      and we can safely call CTRSV if 1/M(n) and 1/G(n) are both greater than
  210.      max(underflow, 1/overflow).
  211.  
  212.  
  213.  
  214.  
  215.  
  216.  
  217.  
  218.  
  219.  
  220.  
  221.  
  222.  
  223.  
  224.  
  225.  
  226.  
  227.  
  228.  
  229.  
  230.  
  231.  
  232.  
  233.  
  234.  
  235.  
  236.  
  237.  
  238.  
  239.  
  240.  
  241.  
  242.  
  243.  
  244.  
  245.  
  246.  
  247.  
  248.  
  249.  
  250.  
  251.  
  252.  
  253.  
  254.  
  255.  
  256.  
  257.  
  258.  
  259.  
  260.  
  261.                                                                         PPPPaaaaggggeeee 4444
  262.  
  263.  
  264.  
  265.